######## REPLICATION SCRIPT FOR "The Influence of Episodic Information on Political Elites: Evidence from Chile" #########
######## Author: Daniel Cruz                                                                          #########

### Packages:
library(tidyverse)
library(stargazer)
library(ggpubr)
library(stats)
library(foreign)
library(haven)
library(labelled)
library(lfe)
library(broom)
library(ggstats)

### Environment
setwd("")

### Load Data
D <- read.csv("Replication Data - Elite.csv")

### FIGURE 1 ####

# a) left panel: T-Test
# _t1 is the treatment with neg. episodic information and positive statistical.
# _t2 is the treatment with pos. episodic informatoin and negative statistical.


DatPol2 <- D %>%
  select(Pol2_t1, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo) %>%
  mutate(Pol2 = Pol2_t1) %>%
  filter(Pol2_t1 %in% c(1:10)) %>%
  bind_rows(D %>% 
              select(Pol2_t2, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo)%>%
              mutate(Pol2=Pol2_t2) %>%
              filter(Pol2_t2 %in% c(1:10))) %>%
  mutate(treat.stat=ifelse(is.na(Pol2_t1) == T, 1, 0))

p2 <- list()
p2[[1]] <- t.test(D$Pol2_t1)
p2[[2]] <- t.test(D$Pol2_t2)

GP2 <- tibble(id=seq(2),
              mean=sapply(p2, function(x) x$estimate),
              conf_lower=sapply(p2, function(x) x$conf.int[1]),
              conf_upper=sapply(p2, function(x) x$conf.int[2])) 

GP2 <- GP2 %>%
  mutate(treatment=ifelse(id==1, "Neg, 1%", "Pos, 20%"))


ggplot(GP2, aes(x=treatment, y=mean))+
  geom_point(color="grey30") +
  geom_pointrange(aes(ymin=conf_lower, ymax=conf_upper), 
                  color="grey30", size=0.5)+
  ggtitle("Policy 2")+
  #ylim(2,7.5)+
  xlab("") +
  ylab("") +
  theme_bw() +
  theme(plot.title = element_text(size=10))


#### Figure 2 ####
# T.test for the second experiment.
# Q6_coded is the outcome variable for the third experiment. Recall that 
# a coder blinded to the treatments was used for this. The accuracy scale ranged 
# from 0 to 2. 
# Q1_Temp==1 idndicate those treated with episodic
# Q1.1_Temp==1 indicate those treated with statistical

t.test(D$Q6_coded[D$Q1_Temp==1], D$Q6_coded[D$Q1.1_Temp==1])

p3 <- list()
p3[[1]] <- t.test(D$Q6_coded[D$Q1_Temp==1])
p3[[2]] <- t.test(D$Q6_coded[D$Q1.1_Temp==1])

G3 <- tibble(id=seq(2),
             mean=sapply(p3, function(x) x$estimate),
             conf_lower=sapply(p3, function(x) x$conf.int[1]),
             conf_upper=sapply(p3, function(x) x$conf.int[2])) 
G3 <- G3 %>%
  mutate(treatment=ifelse(id==1, "Episodic", "Statistic"))



ggplot(G3, aes(x=treatment, y=mean))+
  geom_point(color="grey30") +
  geom_pointrange(aes(ymin=conf_lower, ymax=conf_upper), 
                  color="grey30", size=0.5)+
  ggtitle("Coded Recall")+
  ylim(0.2, 1.2)+
  xlab("") +
  ylab("") +
  theme_bw() + 
  theme(plot.title=element_text(size = 10))





#### Figure 3 ####
# t.test for the third experiment.
# the News_ variable indicate 1 when topic B was chosen and 0 when it wasn't 
# _t1 indicate treatment group with Topic B as Episodic
# _t2 indicate treatment group with Topic B as Thematic

t.test(D$News_t1, D$News_t)

p4 <- list()
p4[[1]] <- t.test(D$News_t1, conf.level = .95)
p4[[2]] <- t.test(D$News_t2, conf.level = .95)

G4 <- tibble(id=seq(2),
             mean=sapply(p4, function(x) x$estimate),
             conf_lower=sapply(p4, function(x) x$conf.int[1]),
             conf_upper=sapply(p4, function(x) x$conf.int[2])) 
G4 <- G4 %>%
  mutate(treatment=ifelse(id==1, "Episodic", "Thematic"))



ggplot(G4, aes(x=treatment, y=mean))+
  geom_point(color="grey30") +
  geom_pointrange(aes(ymin=conf_lower, ymax=conf_upper), 
                  color="grey30", size=0.5)+
  ggtitle("Proportion Favoring Topic B")+
  xlab("") +
  ylab("") +
  ylim(0.3,1) +
  theme_bw() +
  theme(plot.title=element_text(size = 10))


#### Figure 4 ####
## Here I'm comparing elites with citizens. The logic is the same as previously. 
citizens <- read.csv("Replication Data - Citizens.csv")
p_cit <- list()
p_cit[[1]] <- t.test(citizens$Pol2_t1, conf.level = .95)
p_cit[[2]] <- t.test(citizens$Pol2_t2, conf.level = .95)

dtcit_all <- tibble(id=seq(2),
                    mean=sapply(p_cit, function(x) x$estimate),
                    conf_lower=sapply(p_cit, function(x) x$conf.int[1]),
                    conf_upper=sapply(p_cit, function(x) x$conf.int[2])) 

dtcit_all <- dtcit_all %>%
  mutate(treatment=ifelse(id==1, "Neg, 1%", "Pos, 20%"),
         group="Citizens-All")

citizens_educ <- subset(citizens, educ %in% c("Superior universitaria", "Superior no universitaria"))

p_cit_educ <- list()
p_cit_educ[[1]] <- t.test(citizens_educ$Pol2_t1, conf.level = .95)
p_cit_educ[[2]] <- t.test(citizens_educ$Pol2_t2, conf.level = .95)

dtcit_educ <- tibble(id=seq(2),
                     mean=sapply(p_cit_educ, function(x) x$estimate),
                     conf_lower=sapply(p_cit_educ, function(x) x$conf.int[1]),
                     conf_upper=sapply(p_cit_educ, function(x) x$conf.int[2])) 

dtcit_educ <- dtcit_educ %>%
  mutate(treatment=ifelse(id==1, "Neg, 1%", "Pos, 20%"),
         group="Educated Citizens")

dtcit <- rbind(dtcit_all, dtcit_educ)

GP2$group <- "Politicians"

dtcit <- rbind(dtcit, GP2[c("id", "mean", "conf_lower",
                            "conf_upper", "treatment",
                            "group")])

ggplot(dtcit, aes(x=treatment, y=mean, shape=group))+
  geom_pointrange(aes(ymin=conf_lower, ymax=conf_upper), 
                  size=0.5, color="grey30",
                  position=position_dodge(width = 0.6))+
  ggtitle("Influence")+
  ylim(2,7.5)+
  xlab("") +
  ylab("") +
  theme_bw() +
  theme(plot.title = element_text(size=10))


#### Appendix ####
## Per Vignette Data

# 1. Balance Tests

## Balance Exp1
DatPol2 <- D %>%
  select(Pol2_t1, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo) %>%
  mutate(Pol2 = Pol2_t1) %>%
  filter(Pol2_t1 %in% c(1:10)) %>%
  bind_rows(D %>% 
              select(Pol2_t2, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo)%>%
              mutate(Pol2=Pol2_t2) %>%
              filter(Pol2_t2 %in% c(1:10))) %>%
  mutate(treat.stat=ifelse(is.na(Pol2_t1) == T, 1, 0),
         fem_ratio=ifelse(gender=="Femenino", 1, 0),
         unideg_ratio=ifelse(educ=="Superior universitaria", 1, 0),
         local_ratio=ifelse(position %in% c("Alcalde", "Concejal", "CORE"), 1, 0),
         leftw_ratio=ifelse(party_ideo == "Left", 1, 0))


covariates <- c("fem_ratio", "age", "unideg_ratio", "yearsexp_num", "local_ratio", "leftw_ratio")

covariate_means <- DatPol2 %>%
  group_by(treat.stat) %>%
  summarise(female_mean = mean(fem_ratio, na.rm=T),
            female_cilow = t.test(fem_ratio)$conf.int[1],
            female_ciup = t.test(fem_ratio)$conf.int[2],
            age_mean = mean(age, na.rm=T),
            age_cilow = t.test(age)$conf.int[1],
            age_ciup = t.test(age)$conf.int[2],
            degree_mean = mean(unideg_ratio),
            degree_cilow = t.test(unideg_ratio)$conf.int[1],
            degree_ciup = t.test(unideg_ratio)$conf.int[2],
            yexp_mean = mean(yearsexp_num, na.rm=T),
            yexp_cilow = t.test(yearsexp_num)$conf.int[1],
            yexp_ciup = t.test(yearsexp_num)$conf.int[2],
            local_mean = mean(local_ratio),
            local_cilow = t.test(local_ratio)$conf.int[1],
            local_ciup = t.test(local_ratio)$conf.int[2],
            leftw_mean = mean(leftw_ratio),
            leftw_cilow = t.test(leftw_ratio)$conf.int[1],
            leftw_ciup = t.test(leftw_ratio)$conf.int[2])



covariate_means_long <- covariate_means %>%
  pivot_longer(cols = -treat.stat, names_to = c("Covariate",".value"), names_sep = "_") %>%
  mutate(treat.stat = ifelse(treat.stat==0, "Treat 1", "Treat 2"))

ggplot(covariate_means_long, aes(x=mean, y=as.factor(treat.stat), color=as.factor(treat.stat))) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_pointrange(aes(xmin=cilow, xmax=ciup), 
                  , size=0.5) +
  facet_wrap(~Covariate, scales = "free_x") +
  ggtitle("Balance test - Experiment 1")+
  #xlim(-1.2,0.1)+
  ylab("") +
  theme_bw() +
  theme(plot.title = element_text(size=10),
        legend.position = "none") +
  scale_color_grey() 

  ## 2. Balance Exp2
DatPol3 <- D %>%
  select(Q6_coded, Q1_test, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo) %>%
  mutate(fem_ratio=ifelse(gender=="Femenino", 1, 0),
         unideg_ratio=ifelse(educ=="Superior universitaria", 1, 0),
         local_ratio=ifelse(position %in% c("Alcalde", "Concejal", "CORE"), 1, 0),
         leftw_ratio=ifelse(party_ideo == "Left", 1, 0))


covariate_means <- DatPol3 %>%
  filter(Q1_test %in% c(0,1)) %>%
  group_by(Q1_test) %>%
  summarise(female_mean = mean(fem_ratio, na.rm=T),
            female_cilow = t.test(fem_ratio)$conf.int[1],
            female_ciup = t.test(fem_ratio)$conf.int[2],
            age_mean = mean(age, na.rm=T),
            age_cilow = t.test(age)$conf.int[1],
            age_ciup = t.test(age)$conf.int[2],
            degree_mean = mean(unideg_ratio),
            degree_cilow = t.test(unideg_ratio)$conf.int[1],
            degree_ciup = t.test(unideg_ratio)$conf.int[2],
            yexp_mean = mean(yearsexp_num, na.rm=T),
            yexp_cilow = t.test(yearsexp_num)$conf.int[1],
            yexp_ciup = t.test(yearsexp_num)$conf.int[2],
            local_mean = mean(local_ratio),
            local_cilow = t.test(local_ratio)$conf.int[1],
            local_ciup = t.test(local_ratio)$conf.int[2],
            leftw_mean = mean(leftw_ratio),
            leftw_cilow = t.test(leftw_ratio)$conf.int[1],
            leftw_ciup = t.test(leftw_ratio)$conf.int[2])



covariate_means_long <- covariate_means %>%
  pivot_longer(cols = -Q1_test, names_to = c("Covariate",".value"), names_sep = "_") %>%
  mutate(Q1_test = ifelse(Q1_test==0, "Treat 1", "Treat 2"))


ggplot(covariate_means_long, aes(x=mean, y=as.factor(Q1_test), color=as.factor(Q1_test))) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_pointrange(aes(xmin=cilow, xmax=ciup), 
                  , size=0.5) +
  facet_wrap(~Covariate, scales = "free_x") +
  ggtitle("Balance test - Experiment 2")+
  #xlim(-1.2,0.1)+
  ylab("") +
  theme_bw() +
  theme(plot.title = element_text(size=10),
        legend.position = "none") +
  scale_color_grey() 

  ## 3. Balance Exp3

DatPol4 <- D %>%
  select(News_t1, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo) %>%
  mutate(News_exp = News_t1) %>%
  filter(News_t1 %in% c(1:10)) %>%
  bind_rows(D %>% 
              select(News_t2, gender, age, educ, ideo, year_exp, yearsexp_num, position, party_ideo)%>%
              mutate(News_exp=News_t2) %>%
              filter(News_t2 %in% c(1:10))) %>%
  mutate(treat.stat=ifelse(is.na(News_t1) == T, 1, 0),
         fem_ratio=ifelse(gender=="Femenino", 1, 0),
         unideg_ratio=ifelse(educ=="Superior universitaria", 1, 0),
         local_ratio=ifelse(position %in% c("Alcalde", "Concejal", "CORE"), 1, 0),
         leftw_ratio=ifelse(party_ideo == "Left", 1, 0))


covariates <- c("fem_ratio", "age", "unideg_ratio", "yearsexp_num", "local_ratio", "leftw_ratio")

covariate_means <- DatPol4 %>%
  group_by(treat.stat) %>%
  summarise(female_mean = mean(fem_ratio, na.rm=T),
            female_cilow = t.test(fem_ratio)$conf.int[1],
            female_ciup = t.test(fem_ratio)$conf.int[2],
            age_mean = mean(age, na.rm=T),
            age_cilow = t.test(age)$conf.int[1],
            age_ciup = t.test(age)$conf.int[2],
            degree_mean = mean(unideg_ratio),
            degree_cilow = t.test(unideg_ratio)$conf.int[1],
            degree_ciup = t.test(unideg_ratio)$conf.int[2],
            yexp_mean = mean(yearsexp_num, na.rm=T),
            yexp_cilow = t.test(yearsexp_num)$conf.int[1],
            yexp_ciup = t.test(yearsexp_num)$conf.int[2],
            local_mean = mean(local_ratio),
            local_cilow = t.test(local_ratio)$conf.int[1],
            local_ciup = t.test(local_ratio)$conf.int[2],
            leftw_mean = mean(leftw_ratio),
            leftw_cilow = t.test(leftw_ratio)$conf.int[1],
            leftw_ciup = t.test(leftw_ratio)$conf.int[2])



covariate_means_long <- covariate_means %>%
  pivot_longer(cols = -treat.stat, names_to = c("Covariate",".value"), names_sep = "_") %>%
  mutate(treat.stat = ifelse(treat.stat==0, "Treat 1", "Treat 2"))

ggplot(covariate_means_long, aes(x=mean, y=as.factor(treat.stat), color=as.factor(treat.stat))) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_pointrange(aes(xmin=cilow, xmax=ciup), 
                  , size=0.5) +
  facet_wrap(~Covariate, scales = "free_x") +
  ggtitle("Balance test - Experiment 3")+
  #xlim(-1.2,0.1)+
  ylab("") +
  theme_bw() +
  theme(plot.title = element_text(size=10),
        legend.position = "none") +
  scale_color_grey() 

## 2. Full Models


# Models for experiment 1 
model1 <- lm(Pol2 ~ treat.stat + yearsexp_num, data=DatPol2)
model1.1 <- lm(Pol2 ~ treat.stat + yearsexp_num + ideo, data=DatPol2)
model1.2 <- lm(Pol2 ~ treat.stat + yearsexp_num + ideo + educ, data=DatPol2)
model1.3 <- lm(Pol2 ~ treat.stat + yearsexp_num + ideo + educ + gender , data=DatPol2)

stargazer(model1, model1.1, model1.2, model1.3,
          title="Results with Covariates - Exp 1.2",
          type = "text",
          float = TRUE,
          no.space = TRUE,
          header=FALSE,
          font.size = "small",
          intercept.bottom = F,
          digits = 2,
          t.auto = F,
          p.auto = F,
          notes.align = "l",
          notes.append = TRUE,
          omit.stat=c("f", "ser"),
          suppress.errors = FALSE,
          column.sep.width = "2pt"
)

## Interaction Models Exp 1.
DatPol2$Experienced <- ifelse(DatPol2$yearsexp_num > 1, 1, 0)
experience_int <- lm(Pol2 ~ treat.stat*Experienced + ideo + educ + gender , data=DatPol2)
gender_int <- lm(Pol2 ~ treat.stat*gender + yearsexp_num + ideo + educ , data=DatPol2)
ideo_int <- lm(Pol2 ~ treat.stat*party_ideo + yearsexp_num + educ + gender , data=DatPol2)
DatPol2$position_dic <- ifelse(DatPol2$position %in% c("Alcalde", "CORE", "Concejal"), "Local", "National")
position_int <- lm(Pol2 ~ treat.stat*position_dic + yearsexp_num + ideo + educ + gender , data=DatPol2)

stargazer(experience_int, gender_int, ideo_int, position_int,
          title="Interaction Models - Exp 1",
          type = "text",
          dep.var.labels = "Policy Evaluation",
          float = TRUE,
          no.space = TRUE,
          header=FALSE,
          font.size = "small",
          intercept.bottom = F,
          digits = 2,
          t.auto = F,
          p.auto = F,
          notes.align = "l",
          notes.append = TRUE,
          omit.stat=c("f", "ser"),
          suppress.errors = FALSE,
          column.sep.width = "2pt"
)

## Models for experiment 2 
model2 <- lm(Q6_coded ~ Q1_test + yearsexp_num, data=D)
model2.1 <- lm(Q6_coded ~ Q1_test + yearsexp_num + ideo, data=D)
model2.2 <- lm(Q6_coded ~ Q1_test + yearsexp_num + ideo + educ, data=D)
model2.3 <- lm(Q6_coded ~ Q1_test + yearsexp_num + ideo + educ + gender , data=D)

stargazer(model2, model2.1, model2.2, model2.3,
          title="Results with Covariates - Exp 2",
          type = "text",
          float = TRUE,
          no.space = TRUE,
          header=FALSE,
          font.size = "small",
          intercept.bottom = F,
          digits = 2,
          t.auto = F,
          p.auto = F,
          notes.align = "l",
          notes.append = TRUE,
          omit.stat=c("f", "ser"),
          suppress.errors = FALSE,
          column.sep.width = "2pt"
)


## 3. Intercoder Reliability Test
library(tidycomm)

coder1 <- read.csv("q6_coder_1.csv")
coder1 <- coder1 %>%
  mutate(coder = "coder1",
         text.id = 1:nrow(coder1))

coder2 <- read.csv("q6_coder_2.csv")
coder2 <- coder2 %>%
  mutate(coder = "coder2",
         text.id = 1:nrow(coder2))

intercoder <- rbind(coder1, coder2)

int_rel <- intercoder %>% as_tibble %>%
  test_icr(text.id, coder, code, cohens_kappa = T)

xtable::xtable(int_rel)


